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Abstract 

Quadrature-based moment-closure methods are a class of approximations that 
replace high-dimensional kinetic descriptions with lower-dimensional fluid mod- 
els. In this work we investigate some of the properties of a sub-class of these 
methods based on bi-delta, bi-Gaussian, and bi-B-splinc representations. We 
develop a high-order discontinuous Galerkin (DG) scheme to solve the result- 
ing fluid systems. Finally, via this high-order DG scheme and Strang operator 
splitting to handle the collision term, we simulate the fluid-closure models in the 
context of the Vlasov-Poisson-Fokker-Planck system in the high- field limit. We 
demonstrate numerically that the proposed scheme is asymptotic-preserving in 
the high-field limit. 

Keywords: Asymptotic-Preserving; Discontinuous Galerkin; Vlasov-Poisson; 
Fokker-Planck; Moment-Closure; Moment-Realizability; Plasma Physics; 
High- Order Schemes 



1. Introduction 

The focus of this work is on ID fluid models of plasma with a class of fluid- 
closure approximations known as quadrature-based moment- closures. In particu- 
lar, in this work we are interested in applying these fluid-closure approximations 
to the one-dimensional form the Vlasov-Poisson-Fokker-Planck (VPFP) equa- 
tions: 

~fi + i~f.i--Elv = ti{vf + ^~f.,\ , E.,^ — {~po^~p), (1) 
m \ m J - meo 
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where t G R is time, i G K is the spatial coordinate, w € M is the velocity, 
f{t,x,v) is the probability density function for electrons, E(t,x) is the electric 
field, and p ^ J m fdv is the electron mass density. The parameters in this 
equation are the elementary charge, e, the electron mass, m, the Boltzmann 
constant, ks, the temperature of the equilibrium state, 0, and the stationary 
background ion mass density, po{x). In the above expression ~ is used to denote 
dimensional dependent and independent variables (i.e., these decorations will be 
removed after non-dimensionalization) . These equations describe the dynam- 
ics of electrons (as represented by the PDF f{i,x,v)) that evolve via Coulomb 
interactions and collisions in the form of a Fokker-Planck drift-diffusion opera- 
tor. The Fokker-Planck operator tries to drive the system to a thermodynamic 
equilibrium with constant temperature 9. 

1.1. The Vlasov-Poisson-Fokker-Planck system in the high-field limit 

Fluid-closure methods as described in this work will generally not accurately 

approximate solutions of ([T|) in the collisionless limit, £ — oo. Instead, we focus 

here on the high-collision limit; and in particular, the high- field limit, which 

describes the long-time, large-scale, high-coUisional, and large electric field limit 

of the VPFP system. 

In order to derive the high-field limit of the VPFP system we introduce 

a non-dimensionalization via the characteristic scaling constants: T (time), L 

(length), Eq (electric field), and N (number density), such that 

i = Tt, i = Lx, i> = LT-'^v, E = EqE, p = mNp, f NTL^^f. 
This reduces the VPFP system (H]) to 

^- - - m) -K"^ - (^) ^■•) ,. • ^ (s) - ■ 

We define the dimensionless parameter e = (pT) ^ and choose T, L, and Eq as 
follows: 

T-'-^, L^V^{^), Eo = ^^^V^. (2) 

After simplification, this results in the following non-dimensional Vlasov-Poisson 
Fokker-Planck (VPFP) system: 

f,t + vf,. = -(F{F-'f) ) , i?,, =po-p, (3) 

where p — J fdv and i^(t, v) is the isothermal equihbrium distribution: 

F{t,x,v) = ^cxp (^-ljv + Eit,x) f'^ . (4) 

The so-called high-field limit is when e — > 0+ , which, under the choices of the 
characteristic time, length, and electric field chosen in ([2]), describes long time 
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(T — >■ oo), large-scale {L — >• oo), and large electric field {Eo — oo) dynamics of 
the VPFP system. In particular, Nieto et al. [1^ proved that as e — >■ 0+, the 
solution of the VPFP system converges to the equilibrium distribution 
where (to leading order in e) the mass density, p(t,x), and the electric field, 
E{t,x), satisfy the following non-local advection equation: 

p,t - (pE)^^ = 0, E,., ^po-p. (5) 

1.2. Scope of this work 

Recently, Wang and Jin [3| developed an asymptotic-preserving scheme for 
the VPFP equation, where they modified a fully kinetic solver for VPFP so that 
it remains asymptotic preserving in the high- field limit e — > 0^. The approach 
has the nice property that it can be applied for any value of e > 0. The problem 
with the Wang and Jin scheme is that if one is really interested in regimes 
where e is relatively small (i.e., near thermodynamic equilibrium), then their 
approach is computationally expensive (i.e., requires solving a PDE in 2D rather 
than ID). 

The purpose of this work is to consider fiuid models for the VPFP system 
that not only have the ability to capture the equilibrium dynamics of VPFP 
(i.e., equation ([5])), but also accurately model near-equilibrium dynamics. The 
scope of the current work is to twofold: 

1. We describe and investigate properties of two approaches in the quadrature- 
based moment-closure framework as developed in Fox [l^] (bi-delta distri- 
bution functions) and Chalons, Fox, and Massot [HI (bi-Gaussian distribu- 
tion functions) , and describe a modification of these based on bi-B-spline 
distributions. This is described in Sj2j 

2. We then take this class of quadrature-based moment-closure approaches 
and, via a high-order discontinuous Galerkin scheme with Strang operator 
splitting for the collision operator, approximately solve the VPFP system 
in the high-field limit. The numerical method is developed in S|3] and 



applied to two test problems from Wang and Jin [IJ] in 



2. Quadrature-based moment-closure methods 

In this section we describe three variants of the quadrature-based moment- 
closure approach. In fj2T] we describe the simplest of these: quadrature based 
on delta distributions. Delta distributions have been used in several previous 
works including in gas dynamics applications in Fox [lol | and Yan and Fox 
[2^, multivalued solutions of Euler-Poisson and multiphase solutions of 



the semiclassical limit of the Schrodinger equation [ll|, [l3| . In we describe 
a generalization of this approach using Gaussian distributions developed by 
Chalons, Fox, and Massot |5[. Finally in ij2.3l we propose a modification of the 



Gaussian distribution approach that uses compactly supported B-splines. 
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2.1. Quadrature using delta distributions 

In the quadrature-based moment-closure approach using delta distributions 
we assume a PDF that is a sum of delta functions with unknown weights and 
positions (we show here the simplest version of this approach using only two 
quadrature points): 

fit, X, V) = LJlS {V - 111) + L025 [v - H2) , (6) 

where the parameters uii, /ii, a;2, and /i2 are all functions of t and x. This 
approach is reminiscent of other discrete velocity models such as the Broadwell 
model 0,0; however, a key difference is that the discrete velocities, fii and /X2, 
are potentially different at each point in space and time. 
The first four moments of ([5]) are 

Mq = p = wi/i? + a;2/i2, (7) 

Ml = pu = wi/ij + a;2At2, (8) 

M2 = pu^ +p = ujifij + L02IJ-I, (9) 

M3 = pu^ + 3pu + q = ujifil + a;2M2- (10) 

If we assume that p > and p > 0, then the above relationship between the mo- 
ments (Mq, Ml, M2, M3) and the parameters (/ii, ^2, wi, ^2) is one-to-one (see 
discussion below). In the absence of collisions, these moments satisfy equations 
of the form: 

M^,t + Mf+i,, = (11) 

for £ = 0,1, 2, 3. The moment-closure comes from forcing M4 to come from ([6]) 
(rather than letting it be an independent quantity): 

/oo 
V^J dv = UJif4 + UJ2fit- (12) 
-00 

Therefore, the moment-closure problem is reduced to the following: given 
the first four moments of the system, (Mq, Mi, M2, M3), solve system ((7))- p0)) 
to obtain (/ii, /X2, wi, W2), then use these parameters to calculate M4 via (fT2|) . 

Solving system ([7t-(|10|) is equivalent to finding the quadrature points and 
weights for the following weighted Gaussian quadrature rule: 



g{v)w{v)dv « wi .9 (/ii) + W2 ,9 (/i2) (13) 
with weight function w{v) that satisfies 

/oo 
v''w{v)dv^Mk for fc = 0,1,2,3. (14) 
-00 

If we attempt to make this quadrature rule exact with g{v) = 1, v, w^, and v^, 
we arrive at equations ([7 | -(|10 |) . To find the correct Gaussian quadrature rule. 
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we invoke results from classical numerical analysis (e.g., see pages 220-225 of 
Burden and Faires and look for polynomials of degree up to two that are 
orthogonal in the following weighted inner product: 



{g,h)^ := 



l{v) h(v) w{v) dv. 



(15) 



Such polynomials are easily obtained by starting with monomials in v and ap- 
plying Gram-Schmidt orthogonalization with respect to (1151) : 

V'(°)(«) = l, (16) 
ij;'^^\v) ^ V - u, (17) 
^/;(^) (w) = 3ppv^ — {6ppu + 3pq) v + (3ppu^ — + 3upq) . (18) 

The quadrature points pi and p2 are the two distinct real roots of ^(^^(u): 



2p 



2p 



(19) 



Once the quadrature points are known, the corresponding quadrature weights 
can easily be obtained by solving equations ([7]) and (|8]) for the weights: 



Wl, W2 = - ± 



p^q 



2 2v/p2g2 + 4pp3 



(20) 



2.1.1. Flux form of the bi-delta system 

Once pi, p2, wi, and W2 have been computed, we can evaluate the moment 
closure: replace the true M4 with the following: 



-M4 = piu\ + p2 0J2 = pu^ + &py? + '^uq +- — h — 

p p 

From this we can write down the fully closed system: 

U,t + F{U),^ = 0, 

where 

U = [p, pu, pu^ + p, pu^ + 3pu + q] 
is the vector of conserved variables and 



FiU) 
is the flux function. 



pu, pu + p, pu + 3pu + q, pu + 6pu + Auq H 

P 



P_ 

P 



(21) 

(22) 
(23) 

(24) 
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2.1.2. Hyperbolic structure of the bi-delta system 

One can show that system ((22|) with (|23p - (p4)) . assuming that p > and 
p > 0, is a weakly hyperbolic system of PDEs. The eigenvalues of the flux 
Jacobian, F^u, are 

AW=a(2)=^i and A^^) = = (25) 

The two linearly independent eigenvectors are 

rW=r(2) = [l,A.i,M?,A^?]^ and r^'^ = ^ [l, f,,, . (26) 

One can show that the simple waves associated with each of these eigenvectors 
are linearly degenerate: 

Vc/A^^) • rC^) = for fc = 1,2,3,4, (27) 

where Vc/ is the gradient with respect to the conserved variables (|23|) . For a 
detailed analysis of this system see Chalons, Kah, and Massot 01 . 

2.2. Quadrature using Gaussian distributions 

There are two main difficulties with moment-closure based on quadrature via 
delta distributions: (1) the resulting system is only weakly hyperbolic, which 
means that delta shocks P, 0, Q are gcncrically present in the system, and (2) 
a large number of delta functions may be required to g et g ood agreement with 
smooth distributions with large support. Yuan and Fox |25| developed an adap- 
tive Gaussian quadrature strategy to help overcome this problem. A possible 
alternative to using large number of quadrature points was developed Chalons 
et al. who introduced a quadrature moment-closure based on replacing the 
Dirac delta functions with Gaussian distribution functions. In particular, in 
order to simplify the moment inversion equations, each Gaussian is assumed to 
have the same standard deviation. In the case of two Gaussian distributions 
this results in an assumed distribution function of the form: 

^ 71^ [ — ^ J + 71^ [ — ^ ) ' ^''^ 

where cr is the width. In the limit as (t — > we recover the Dirac delta distribu- 
tion moment-closure method. In this section we describe the moment-inversion 
algorithm needed to convert between moments and the parameters in repre- 
sentation ((28|) . Furthermore, we briefly discuss the hyperbolic structure of the 
resulting evolution equations for the moments of ((28)) . 

There are now five free parameters: (^1 , ^2 , 1^1 , W2 , cr) , which can be obtained 
by solving the following moment-inversion problem on the first five moments: 



Mo 


= wi/i? - 


^ W2/i2, 






(29) 


Ml 


= wi/i} - 


Fa;2/i2, 






(30) 


M2 


= wi/ii - 


F W2M2 ^ 




- ^^2/^2) : 


(31) 


M3 


= wi/i? - 


F a;2M2 H 


h 3(7 (wi/ij 


+ ^^2^2) . 


(32) 


M4 


= wi/ii - 


F W2M2 ^ 


h 6cr [uji^l 


+ W2M2) + ('^iMi + '^2^2) 


■ (33) 
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We write this system in terms of primitive variables: 




+ 



3p2(a2 _ 1) 



(34) 
(35) 
(36) 
(37) 



(38) 



P 



where we have introduced the parameter a: 




(39) 



We note tliat the maximum ahowable value of a is clearly 1 (otherwise a would 
be negative). In fact if a = 1, then the bi-Gaussian and the bi-delta func- 
tion representations are equivalent. What is perhaps less obvious is that the 
minimum allowed value of a is zero. In fact, as a — !• 0, the bi-Gaussian repre- 
sentation collapses into a single Gaussian distribution with width given by the 
temperature p/p. Physically reasonable conditions on the primitive variables 
guarantee that the a that comes from solving the moment equations 
above satisfies a G [0,1] (see Theorem 12. ip . The exact formula for computing 
the value of a is given in Algorithm [TJ 

Theorem 2.1 (Moment-realizability condition). Assume that the primitive 
variables satisfy the following conditions: 

• Positive density: < p, 

• Positive pressure: < p, 

3 I 2 

• Lower bound on r: ^ '"^ < r, 



1. If q =/= then there exists a unique a G (0,1] that satisfies the following 
cubic polynomial: 



Furthermore, from this a we can uniquely obtain the quadrature abscissas 
and weights in order to exactly solve system p4p - psp ; 



• If q = 0, bound on r: 



2 r, 2 



p — — p 



V{a) = 2p^a^ + {pr - 3/) pa - pq^ = 0. 



(40) 




(41) 



(42) 
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2 Q 

2. If q = and ^ < r < then there exists a unique a E (0, 1] such that 



'fc^. (43, 



Furthermore, in this case the quadrature abscissas (j4ip and weights (I42p 
are again the unique solutions of system (j34p - (|38p . 

If q = and r = i/ien a = 0. T/iis case corresponds to a single 
Gaussian distribution. In this case we lose uniqueness of the quadrature 
abscissas and weights, but without loss of generality we can take 

p,i,p,2^u, (44) 
t^i, ^^2 = ^, (45) 



and still exactly solve system ([34]) -(j38]) . 
Proof. We take each point in turn. 

1. Let us momentarily assume that a is known and that a G (0, 1]. In this 
situation we are left with four unknowns: (/ii, /i2, wi, a;2), which are deter- 
mined by satisfying the first four moment equations: ([M|) - ((57)) . Just as in 
the case of quadrature-based moment-closures using delta distributions, 
these equations can be solved by constructing a set of polynomials that 
are mutually orthogonal in the inner product p5|). The weight function 
satisfies: 

if fc == 0, 
q ii k — 3. 

Up to degree 2 these mutually orthogonal polynomials can be written as 

V'(°'(«) = l, (47) 
iP'^^'> (v) = V ~ u, (48) 




i,C^){v) =v' -\2u+^\v+{u' + ^-^\. (49) 
V pa J \ pa p J 

It is easy to show that the two real roots of -0(^)(i;) are (HH). 

The weights (|42]) are obtained by plugging (gl]) into (jMl) and ([35]) and 
solving the resulting 2x2 linear system for wi and lo2- 

Finally, we must obtain a formula for a e (0,1]. This is achieved by 
plugging (PT|) and (|^^ into the final moment equation: After sim- 

plification this yields the cubic polynomial equation given by (PO]). We 
note that under our assumptions we find that 

V{0) = -pq^ < and V{1) = pp(r- ^!_L^^ > o. 

V PP J 
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Therefore by continuity of V{a) we are guaranteed that there exists at 
least one root in (0, 1]. To estabhsh that there is a unique root in (0, 1] 
we note that V{a) is convex in (0, 1]: 

V"{a) ^12p^a>0 in (0,1]. 

The exact formula for computing the value of a is given in Algorithm [T] 

2. li q = and ^ < r < ^ we note that (gO]) reduces to 

P{a)=a {2p^a^ + {pr - 3/)) = 0. 

The unique root of V{a) in (0, 1] is given by With this value of a 

we can again solve (|34p - ([55|) using the /ii, ^2, wi, and u!2 given by (|4ip 
and dUl). 

3. If g = and r = ^y- the only solution of (j40]) is a = 0. In this case the 
moment equations (|34 l) - (p8| reduce to 

Affe = wi^i + for fc = 0,1,2,3,4. 

This system has an infinite number of solutions of the form: 

A'l = = u and cji + CJ2 = P- 

Without loss of generality we take (l44l) and p5|) . 

□ 

2.2.1. Flux form of the hi- Gaussian system 

After obtaining a, /ii, /i2, wi, and 0^2, we impose the following moment- 
closure: 

Afs = /ii + fi2 <A + lOcrMs - Xha^Ah = pu^ + lOpu^ H ^- h aMj, (50) 

where 

Af5 = f 4 + ^ + 10,-.^ + 1M)_ ^ (4,~+ 5pu) , (51) 
\P P P J P 

q = q/a (we set (7 = if a = 0). (52) 

Finally, we write the closed bi-Gaussian quadrature-based moment-closure sys- 
tem in the form (|22p . where 

U ~ [p, pu, pu^ + p, pu'^ + ipu + q, pu^ + 6pu^ + Aqu + , (53) 
F{U) = [pu, pu^ + p, pu^ + 3pu + q, pu"^ + &pu^ + Aqu + r, M5] ^ . (54) 
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Algorithm 1 Formula for computing a for the bi-Gaussian closure model 
1: input: {p, p, q, r) . 

2: 

check: p > 0, p > 0, r> ^+2^. 



if q = then 



if r < ^ then 



9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 



return a = J^l^. 



// algorithm is finished 



else 

return a = 0. // algorithm is finished 
end if 



end if 
let: Q 



o _ pr 



7^ 



£2_ 

4p3 



2? = 7^2 



if 2? > then 

T* = n-Vv. 

if r* < then 

else 
end if 



return a = \/7?. + VP + T. // algorithm is finished 

else 

e = icos-i (7^/yQ3). 

ai = 2VQcos(6I). 
if < ai < 1 then 

return a = ai. // algorithm is finished 

else 

02 = 2x/Qcos(6' + 27r/3). 
if < 02 < 1 then 

return a = 02. // algorithm is finished 

else 

02 = 2v^cos(6l + 47r/3). 

return a = 03. // algorithm is finished 
end if 
end if 
end if 
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2.2.2. Hyperbolic structure of the hi-Gaussian system 

The flux Jacobian of the bi-Gaussian system described above can be written 
in the following form: 






1 

















1 

















1 

















1 


5, Mo 


M5,Mi 






-^'-^53/4 



where U = (Afo, ^'A, M2, A/3, ^^4)^ and A/5 is given by ([5 
and right eigenvectors of the flux Jacobian are of the form 

A^^-) = zu and r'^) = [l, Zk. zl zl zt] 
and the left eigenvectors can be written as 



(55) 

The eigenvalues 
(56) 



5 3,4,5 4,5 5 

i=j+i e=j+i 

m=l+l 



(57) 



n - zj) 



where all of the above sums and products exclude the index k. Therefore, the 
eigenstructure hinges on the values of the five numbers: Zk for k = 1,2,3,4,5. 
Unfortunately, we are not able obtain these quantities in closed form. However, 
it is possible to calculate approximate values for these quantities in certain 
asymptotic limits. In particular, one important limit that is useful later on in 
this work is the near thermodynamic equilibrium limit, which is characterized 
here by fii ~ fi2- In this limit the five distinct z^'s are given by 



21,2:2 



Z3 



24,25 



2q 
5p 
2q 
5p 

2q 
5p 



(5±. 


/To) 39(1 -a) 




P 




-Mil^), 




/lO)p(l-a) 


P 



0(lM2-/ilP) 



0(lM2-/il| 



(58) 
(59) 

(60) 



Assuming that p > and p > and noting that a w if /ii ~ ^2, we see 
that in this limit the system is strongly hyperbolic. We can also approximately 
compute 



CkQP {I - a)^(14p'^a'' + 5pq 
20p(4p3a3 + pq'^y 



0(IM2-Mi| 



(61) 



11 



where 

cfc = |2 + VTo, 2 - VTo, ^,2- VTo, 2 + VTo 

We note that in this limit VfyA''^'^ • A''^ changes sign if q changes sign, which 
means that there must exist some state, U*, for which V^A''""-' • A''^ — 0. This 
shows that the waves in this system are non- genuinely nonlinear and may admit 
composite wave solutions (e.g., see [Tsj). 

In order to illustrate how the bi-Gaussian behaves, we show numerical so- 
lutions of two Riemann problems using the discontinuous Galerkin scheme as 
described in ^ 33. 21 - ^3. 31 . In Figure[l]we show the solution to a Riemann problem 
with a heat flux, q, that is striclty positive. In this problem the solution is a 
set of five classical waves (i.e., shocks and rarefactions). In contrast, in Figure 
[2] we show the solution to a Riemann problem with heat flux, q, that changes 
sign over the simulation domain. In this case we see a solution with two com- 
pound waves, each of which is a rarefaction connected to a shock, propagating 
in opposite directions. 

2. 3. Quadrature using C'' B-splines 

From the above discussion we can view the bi-delta and the bi-Gaussian 
closure methods as members of the same family of methods, where the bi-dclta 
is one extreme (compactly supported distributions with zero width) and the 
bi-Gaussian is the other (non-compactly supported distributions with maximal 
standard deviation as allowed in this representation). Using this point-of-view, 
we can also construct methods that are in-between these two extremes. One 
simple example of this is a bi-distribution representation based on B-splines: 

7(t, X, v) = u:iBl {v - fii) + UJ2BI {v - /.2) , (62) 

f {2v + V^) \l - V^ <2v <0, 
B»^{l{V^^'^v) iiO<2v<V^, (63) 



where 



otherwise. 



If we introduce the parameter 



= ^ (1 - a) , (64) 
P 

and force ((62)) to have as its first five moments Mq, Mi, AI2; M^, M4, then we 
again arrive at system p4|) -p7 |) . but now with a slightly modified equation to 
match the M4 moment: 

UJijii + UJ2P-) = pu + iqu + 6apu'^ + r ^ (3a + 2)(q; — 1). (65) 

5p 

Theorem 2.2 (Moment-realizability condition). Assume that the primitive 
variables satisfy the following conditions: 
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Density att - 0.15 



Pressure at t = 0.1 5 



(a) 



1 - 

0.5- 




(e)- 



0.25 0.5 



Heat flux att-0.15 




Abscissas at t ^ 0.1 5 



(b) 



(d)"^0 



(f) -0^ 



0.25 0.5 



0.25 0.5 0.75 

Weights att -0.15 



Figure 1: A shock tube problem for the bi-Gaussian system. In this example the initial states 
are (p, it,p, 5, r)ij,ft = (1.5,-0.5,1.5,1.0,4.5) and (p, m, p, q, ?')right = (1.0,-0.5,1.0,0.5,3.0). 
This data is chosen so that q > OVx,i, ensuring that we do not encounter points where the 
convexity changes. Shown in these panels arc the (a) density (p), (b) pressure (p), (c) heat flux 
{q), (d) width parameter (o), (e) quadrature abscissas {/ii, (12), and (f) quadrature weights 
(ojl, L02). The resulting solution shows, counting waves from left to right, a 1-rarefaction, 
2-shock, 3-shock. 4-rarcfaction, and 5-shock. 
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Density att - 0.12 Pressure at t - 0.12 

i i 1 3| i i 




Figure 2: A shock tube problem for the bi-Gaussian system. In this example the initial 
states are (p, u, p, q, r)ieft = (1,1,1/3,0,0.3) and (p, fvp, <?, r')right = (1,-1,1/3,0,0.3). This 
problem is similar to the one shown in Chalons et al. Q; however, we have reduced the initial 
value of r from 1/3 to 0.3 in order to emphasis the 1-shock and 5-shock. Shown in these 
panels are the (a) density (p), (b) pressure (p), (c) heat flux (q), (d) width parameter (a), (e) 
quadrature abscissas {ni, P2), and (f) quadrature weights (cji, id2)- The resulting solution 
shows a 1-shock connected to a 2-rarefaction, as well as a 4-rarefaction connected to a 5-shock. 
We note that having compound waves is typical of systems with non-convex fluxes (see for 
example pages 350—357 of LeVeque [T^ . 
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• Positive density: < p, 



• Positive pressure: < p, 

• Bound on r: 3l + sl < r < 1^ + ^ . 

p p — — 3p 13p 

Then there exists a unique a G l] ^^0,^ satisfies the following cubic polyno- 
mial: 

V{a) = IS/a^ - Qp^a^ + ap {brp - I2p^) - bpq^ = 0. (66) 

Using this a and the definitions for the abscissas and weights given by (j41[) and 
we can exactly solve the moment inversion problem given by equations 
dMD-dSTD and dssi). 

Proof. The abscissas and weights given by (jH]) and automatically satisfy 
((M|) - ([57)) for any a G (0, 1]. Therefore, the only thing left to do is to satisfy 
equation (p5)) . Using (PT|) and (|^^ . one finds that (p5)) reduces to the cubic 
polynomial (|66| . Wc compute the following: 



and V"{a) > Va > ^, 

which completes the proof. □ 

Remark. Because the bi-B-spline ansatz is compactly supported, there is a max- 
imum value of r that can represented. In particular, if Mq through M4 represent 
the moments of a Gaussian distribution, the value of r = ip^ / p will exceed the 
maximum allowed value in the above theorem. In practice we can remedy this 
situation by taking a = (i-e., the minimum allowed value) whenever r exceeds 

the maximum allowed value of ^ — h 

3p 13p 

The moment closure imposed by the bi-B-spline ansatz is the following re- 
placement for the Ms moment: 

M5 = pu^ + lOu^ [pu + q) + — {5- 4a) 

2 2 3 i^'^) 

P u , „ , „ 9\ bq u q-^ 
4- ^ (12 + 6a - 13a2) + + 

p pa p-'a-' 

The bi-B-spline moment-closure has two advantages over the bi-Gaussian sys- 
tem: (1) a is uniformly bounded away from 0, which means we don't have to 
worry about the same degeneracies as in the bi-Gaussian case; and (2) the dis- 
tribution function f{t,x,v) is compactly supported and piecewise linear, which 
makes computing integrals such as those needed in the flux-vector splitting 
method described in i)3. 31 simpler. In Figure |3] we show a simulation using the 
bi-B-spline moment-closure on the same initial data as used in Figure [T] 
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Density at t - 0.15 Pressure at t - 0.15 




Figure 3: A shock tube problem for the bi-B-sphnc system. In this example the initial states 
are (p, tt,p, g, r)ioft = (1.5,-0.5,1.5,1.0,4.5) and (p, m, p, q, r')right = (1.0,-0.5,1.0,0.5,3.0). 
Shown in these panels are the (a) density (p), (b) pressure (p), (c) heat flux (g), (d) width 
parameter (a), (e) quadrature abscissas (pi, p2), and (f) quadrature weights (oji, 102)- The 
resulting solution shows, counting waves from left to right, a 1-rarefaction, 2-shock, 3-shock, 
4-rarefaction, and 5-shock. 
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3. DG quadrature-based moment-closure schemes for VPFP 



Wc describe in this seetion an application of the quadrature-based moment 
approach as described in previous sections to a particular set of equations 
from plasma physics: the Vlasov-Poisson-Fokker-Planck (VPFP) equations (|3]). 



Wang and Jin IJ] recently deyeloped an asymptotic-preserving scheme for the 
VPFP equation, where they modified a fully kinetic solver for VPFP so that 
it remains asymptotic preserving in the high- field limit e — >■ 0^. Although the 
Wang and Jin approach has the nice property that it can be applied for any 
value of e > 0, if one is really interested in regimes where s is relatively small 
(i.e., near thermodynamic equilibrium), then their approach is computationally 
expensive (i.e., requires solving a PDF in 2D rather than ID). Our focus in this 
section is on approximately solving the VPFP system using quadrature-based 
moment-closure techniques that remain asymptotic-preserving in the high-field 
limit e — > 0^. This approach allows us to efficiently capture near thermody- 
namic equilibrium solutions. 

3.1. Strang operator splitting 

Wang and Jin [l^ achieve an asymptotic-preserving scheme through the 
use of a clever semi-implicit time discretization. In this work we make use of 
a more standard trick: Strang operator splitting which has been used for 
Vlasov-Poisson simulations since the work of Cheng and Knorr . In particular. 



Schaeffer [23[ modified the Cheng and Knorr approach to construct an efficient 
method for VPFP. 

In our approach we use a Strang splitting for the VPFP system under a 
quadrature-based moment-closure with the following steps: 
1. Solve the Poisson equation: 

-(f>,x.x = Poix) ~ p'^ix), £'" = -0_^. 



2. On [t", t" + ^] and for £ = 0, 1, 2, 3, 4 solve VPFP with only the collision 
operator (Af" ^ M"): 

MLt = ^ {({i - l)Me^2 - lE'^Mi^i ~ IMt) . 

3. On [t", -f At] solve the collisionless quadrature-based moment-closure 
system ([^ with the appropriate definitions for U and F{U) (M" — >■ M"+^). 

4. Solve the Poisson equation: 

-^,.,. = Po{x)~p''+Hx), E"+^ 
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5. On [i" + ^, i" + At] and for £ = 0,1, 2, 3, 4 solve VPFP with only the 
collision operator ^ Af"+^): 



Mi, 



i |£(£ - l)Af£_2 - eE"+^Mi_i - me} 



The spatial discretizations are handled via a high-order discontinuous Galerkin 
discretization, which is briefly described below. 

3.2. High-order discontinuous Galerkin spatial discretization 

We make use of the discontinuous Galerkin (DG) method as developed by 

Cockburn and Shu 0] and implemented in the DoGPack software package [IJ] 

to solve hyperbolic conservations of the form (P^ . 

We begin by constructing an equally spaced numerical grid on [a, b] consisting 

of M elements, each of the form: Ti = [xi — ^,Xi + -^j , where Ax = {h — 

a)/M is the grid spacing. Note that Xi denotes the center of element %. Next 

we define the broken finite element space 

meaning that on each element 7i, v^'^ will be a polynomial of degree at most 
k. The solution, U^'^ G V"^^, restricted to element Ti can be written as 

A; 

U^\~Y.^\t)M£,), (69) 

1=1 

where on each element x = Xi-\-^ (Aa;/2), and (/'(O s-'^s the orthonormal Legendre 
polynomials: 

^(C) = |l,x/3e, ^(3e^-l),...|. (70) 

In order to obtain the semi-discrete DG method we multiply (|22p by fj{C), 
integrate over a single element %, replace the exact U by ((69|) . and integrate- 
by-parts in X. After simplification, this results in the following set of coupled 
ordinary differential equations in time: 

jt^''=h f/iu)^Md^-^[p,im+i-vji-m-^], (71) 

where J-,_ i is the numerical flux at interface x ~ .t,- i , which must be calcu- 

2 2 

lated from an (approximate) Riemann solver (see ^3.3l below'). In order to time 
advance this semi-discrete scheme, we make use of the standard third-order to- 
tal variation diminishing Runge-Kutta (TVD-RK) as described in Gottlieb and 



Shu [12| . We make use of the moment- limiters described in Krividonova [15| to 
suppress unphysical oscillations when required. 

Finally we note that the Poisson equation in the operator split scheme de- 
scribed above is solved using a local discontinuous Galerkin scheme that is 



described in detail in Rossmanith and Seal [22 1. 
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3.3. Kinetic-based Riemann solvers 

One missing ingredient from the discussion of the discontinuous Galerkin 
scheme in the previous section is a description of how the numerical flux, i , 
is computed. Since we have the ability to reconstruct the distribution function 
u) for any we can use a kinetic flux-vector splitting approach (see 

for example Mandal and Deshpande [l^). In kinetic flux- vector splitting we 
split the flux into right-going contributions immediately to the left of interface 

1 and left-going contributions immediately to the right of interface x^_^^. 



^r-h=^tl+^^L^ (72) 



where 



= j f(t,x._i,vj dv and ^i_i= J dv. (73) 

3.4- Stiff source term solution and the asymptotic-preserving condition 

The final missing part of the Strang operator split algorithm presented in 
ij3.1l is the solver for the collision operator. A big advantage of considering 
fluid solvers over kinetic solvers in the context of VPFP is that the diffusion 
operator in v becomes an ODE for the moments. In particular, in the Strang 
split approach detailed in §3.11 the electric field, E{t, x). is frozen in time during 
each of the collision operator steps, meaning that the resulting ODEs are linear 
constant coefhcient equations that can easily be solved exactly. The full solution 
over a time step [t", t" -f M] with initial data (Mg", Mf , M^, M3", M4") is 

Af|^'+i = Mg", (74) 
Mf +1 = Z {EAq + Afi") - EM^, (75) 



2EZ {EM^ + + (1 + E^) M^ 



' 



M3"+i = Z^ {E {E^ - 3) A/(5' + 3{E^- l) + 3EM^ + MJ}) 

- 3EZ^ {{E^ - 1) M'a' + 2EMI' + M^) (77) 
+ 'i{E^ + \)Z {EM^ + Ml') -E{E^ + 3) 

M2+^ = Z^{{E'^ - 6E^ + 3) A/^' + AE {E^ - 3) A/f -f {E^ - l) 
+ AEAl^' + AI2) - AEZ^ {E {E^ - 3) A/^' + 3AII'{E^ - l) 
+ 3EAI2 + AI^) +6{E^ + 1) Z^ {{E^ - l) AI^ + 2EAI{' + A'^ 

- AE {E^ + 3) Z (EAl^' + M[') + {E"^ + 6E^ + 3) Ad^, 

where £; = and Z cxp [-At/e]. 



(78) 
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4. Numerical simulations in the high-field limit 

In order to verify the proposed DG operator split method using the quadrature- 
based moment-closure we consider two test cases from Wang and Jin [31 : (1) 
verification of the asymptotic-preserving property and (2) a periodic Riemann 
problem. All of these problems are defined on [0, 1] with periodic boundary 
conditions. 



Verification of the asymptotic-preserving property 
In order to verify the asymptotic-preserving property of the proposed scheme, 
we attempt two versions of the same problem from Wang and Jin 14 [ . 
For the first problem we start with an isothermal Gaussian: 



f{0,x,v) 



2tt 
2n 



exp 



4<» 



(2-f cos {2ttx)). 



with a neutralizing background charge of 

^2^ 



Po{x) 



; cxp [cos(27ra;)] 



(79) 
(80) 

(81) 



1.2660658777520083 
The quantity \\Mi+ pE\\i2 is plotted for various e as a function of time in Figure 
mja). These resuhs verify that ||Mi + pEW^i = O (e) for aU t. 

For the second problem we start with initial data that is not in isothermal 
equilibrium: 



f{0,x,v) 



p(0,x) 



2V27r 



cxp 



1 



(f + 1.5)' 



cxp 



1 



(82) 



with the same p{Q^x) and neutralizing background charge as in the previous 
example. The quantity || Mi -I- pi? 11^,2 is plotted for various e as a function of time 
in Figure HJb). These results verify that the numerical schemes immediately 
drives the non-equilibrium initial data near the equilibrium distribution such 
that \\Mi + pE\\l2 =0{e). 

The simulations presented in Figure|4]were done with the bi-Gaussian moment- 
closure. The bi-delta and bi-B-spline methods give near identical results for this 
problem. 

4.-2. Double periodic Riemann problem 

The initial data is the distribution function (|79p with 

f (1/8, 1/2) ifO<a;<l/4, 
(p(0, x), po{x)) = \ (1/2, 1/8) if 1/4 < x < 3/4, (83) 
[(1/8,1/2) if3/4<x<l. 

The solution using the bi-B-spline moment-closure is shown in Figure [51 The 
results agree well with those in Wang and Jin IJ]. The bi-Gaussian moment- 
closure has difficulties with this problem due to the steep gradients in the solu- 
tion in regions where a is small but non-zero. 
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||M^+p E||^ vs. time 



||M^+p E||^ vs. time 



-eps=1e-3 
-eps=1e-4 
-eps=1e-5 
eps=1e-6 



^QQQOOOOOQOOOOOOOOOOOOOOOOOOOQOOOOOOOOOOOOOOOO 



(a) 



xmxoocoaxoocxxxiC i 



-©-eps^ 


1e 


-3 


-H-eps= 


1e 


-4 


-0-eps^ 


1e 


-5 


— M— eps= 


1e 


-6 



(b) 0^ 



Figure 4: Verification of the asymptotic preserving capability of the method. Shown in these 
two plots are the L^-norms of the difference in the fiuid mean velocity (u) and the equilibrium 
mean velocity {—pE) as a function of time. In Panel (a) the initial conditions arc already 
in equilibrium, while in Panel (b) the initial conditions is not in equilibrium. We show only 
every 20th time step value so that the various curves are more easily identified. We note that 
for e = 10~^ and e = 10^^, the resolution in the simulations is not enough to resolve the 
non-equilibrium deviations; and therefore, the difference between u and —pE is negligible. All 
runs were done with 64 mesh elements. 



Density at t- 0.2 M|att = 0.2 




Figure 5: Periodic Riemann problem of Wang and Jin Il4ll . Shown arc the solutions with 
100 elements (blue circles) and 2000 elements (red line). This simulation is difficult for the 
bi-Gaussian moment-closure due to the steep gradients in the solution in regions where a is 
small but non-zero. These simulations were instead done with the bi-B-spline moment-closiure 
and the results show good agreement with the results of Wang and Jin [l4j . 
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5. Summary 



In this work wc considered quadrature-based moment-closure methods using 
two quadrature points. We briefly investigated the properties of these methods 
and showed connections between bi-delta, bi-Gaussian, and bi-B-spline quadra- 
ture methods. We then applied this formulation to the Vlasov-Poisson-Fokker- 
Planck system in the high-field limit, and, using a high-order discontinuous 
Galerkin scheme with Strang operator splitting, verified the scheme on two test 
problems. Future work will focus on multidimensional plasma physics applica- 
tions. 

Acknowledgements. This work was supported in part by NSF grant DMS- 
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